(** Common histogram implementation *)
open Glob

type info={name:string;title:string}
type geom={n:int;x1:float;x2:float}

module Histogram_axis = struct
  type t = {n_bins: int; max: float; min: float;
  interval: float; rev_int: float; title: string;}

  let print a =
      Printf.printf "n=%i (%f,%f) -%f- /%f : " 
      a.n_bins a.max a.min a.interval a.rev_int;
      print_endline a.title

  let make xa name =
      let interval1 = (xa.x2-.xa.x1)/.float_of_int(xa.n) in
      { n_bins = xa.n; max = xa.x2; min = xa.x1; 
      interval = interval1; rev_int = 1. /. interval1; title = name }
      (*in print a; a*)

  let subst_n a n = make {n=n;x1=a.min;x2=a.max} a.title

  let to_geom a = {n=a.n_bins;x1=a.min;x2=a.max} 

  let name a = a.title
  let range a = a.max -. a.min
  let bounds a = a.min, a.max
  let interval a k = a.interval
  let min a = a.min
  let max a = a.max

  (** centre of bin *)
 let c_of_bin a k =
    a.min +. a.interval *. (0.5 +. float k)
  (** left edge of bin *)
 let l_of_bin a k =
    a.min +. a.interval *. (float k)
  (** right edge of bin *)
 let r_of_bin a k =
    a.min +. a.interval *. (1.+.float k)
  (** transformed coordinate of bin *)
 let w_of_bin a kx =
    a.min +. a.interval *. kx

 let bin_of_c a c = int_of_float (
     (*(x -. a.min) /. a.interval)*)
     (c -. a.min) *. a.rev_int)

 (* TODO : what is "except_last" ??? *)
 let sub a ?(except_last=false) (c1,c2) =
    let n_1_bin = Pervasives.max 0 (bin_of_c a c1) in
    let n_last_bin'' = bin_of_c a c2 - (if except_last then 1 else 0) in
    let n_last_bin' = Pervasives.min n_last_bin'' (a.n_bins-1) in
    let n_last_bin = if n_last_bin' < n_1_bin then n_1_bin else n_last_bin' in
    let n_bins = n_last_bin - n_1_bin + 1 in
    make {n=n_bins; 
       x1 = l_of_bin a n_1_bin; x2 = r_of_bin a n_last_bin} (name a)

  let dots a =
     Array.init a.n_bins (c_of_bin a) 

  let left_dots a =
     Array.init a.n_bins (l_of_bin a) 

 let in_range a x = 
     x >= a.min && x < a.max

  module Bin = struct
    (** CONTAINS SIGN *)
    let ( &: ) = in_range
    let ( |> ) x a = 
      x >= a.max
    let ( <| ) x a =
      x < a.min
  end

 let sgn a c = 
     let range = a.max -. a.min in
     let r = int_of_float ((c -. a.min) /. range) in
     if r>=0 && r<1
     then 1
     else if r>=1 
     then 2
     else if r<0
     then 0
     else failwith "unexpected nan"
end
module Axis = Histogram_axis

exception Incorrect_pair of string

(** functor to make binary operation (plus, divide, etc) 
    for appropriate histogram: H1.hist or H2.hist *)
module Make_bin_op (HS: sig
    type hist
    val get_info: hist -> info (* get info *)
    val dup: ?errors:bool -> info -> hist -> hist (* duplicate hist. with given info *)
    val fit: hist -> hist -> bool (* do axes coincide? *)
    module Array: sig
      type t (* for instance, float array or float array array *)
      type place_info (* for instance, int or int*int *)
      (* note we use syntax a.(p) everywhere: when p has type int or int*int *)
      val get: t -> place_info -> float
      val set: t -> place_info -> float -> unit
      val iterk: (place_info -> float -> unit) -> t -> unit
    end
    val get_data: hist -> Array.t
    (* get squares of errors, if available: *)
    val get_errors: hist -> Array.t option
end) = 
struct
  open HS

  module Aux = struct
    let make_info op info1 info2 =
        { name = info1.name^op^info2.name;
          title = info1.title^" "^op^" "^info2.title; }
    let simple_info s = {name = s; title = s}
    let may errors f =
        ( match errors with
          | None -> ()
          | Some e -> f e
        )
    let may_errors h f =
        may (get_errors h) f
    let is_errors h = 
        match get_errors h with Some _ -> true | _ -> false
  end
  open Aux

  let mul_const h c =
      let info = make_info "*" (get_info h) 
        (simple_info (string_of_float c)) in
      let h' = dup ~errors:(is_errors h) info h in
      let data = get_data h' in
      Array.iterk (fun p v -> 
          data.(p) <- v *. c;
          may_errors h (fun e -> 
              may_errors h' (fun e' -> 
                  e'.(p) <- e.(p) *. c**2.
          ))
      ) (get_data h);
      h'
  ;;
  
(*  exception Bin_op_error of string
  let fail s = raise (Bin_op_error s)*)
  let fail s = raise (Incorrect_pair s)

  let make_binop ~new_info ~calc_data ~calc_error h1 h2 =
      if not (fit h1 h2) then fail "plus: different axes";
      let info = new_info (get_info h1) (get_info h2) in
      let h' = dup ~errors:(is_errors h1 && is_errors h2) info h1 in
      let data = get_data h' and
      data1 = get_data h1 and data2 = get_data h2 in
      Array.iterk (fun p data1_p ->
          data.(p) <- calc_data data1_p data2.(p);
          may_errors h1 (fun e1 ->
              may_errors h2 (fun e2 ->
                  may_errors h' (fun e' ->
                      e'.(p) <- calc_error
                      ~data1:data1_p ~data2:data2.(p) e1.(p) e2.(p)
                  )
              )
          )
      ) data1;

      h'

  let plus ?(c1=1.) ?(c2=1.) h1 h2 = make_binop h1 h2 ~new_info:(make_info "+")
        ~calc_data:(fun data1 data2 -> c1*.data1+.c2*.data2)
        ~calc_error:(fun ~data1 ~data2 sqerr1 sqerr2 ->
            c1**2. *. sqerr1 +. c2**2. *. sqerr2)
  ;;

(*  let mul *)


(*
  let plus h1 h2 =
      if not (fit h1 h2) then fail "plus: different axes";
      let info = make_info "+" (get_info h1) (get_info h2) in
      let h' = dup ~errors:(is_errors h1 && is_errors h2) info h1 in
*)

  let mul h1 h2 =
      if not (fit h1 h2) then fail "mul: different axes";
      let info = make_info "/" (get_info h1) (get_info h2) in
      let h' = dup ~errors:(is_errors h1 && is_errors h2) info h1 in
      let data = get_data h' and 
      data1 = get_data h1 and data2 = get_data h2 in
      Array.iterk (fun p data1_p ->
          data.(p) <- data1_p *. data2.(p);
          may_errors h1 (fun e1 ->
              may_errors h2 (fun e2 ->
                  may_errors h' (fun e' ->
                      e'.(p) <- e1.(p) *. data2.(p)**2. +. 
                                e2.(p) *. data1_p**2.
          )))
      ) data1;
      h'
  ;;

  (* marking by '--' means copypasted from ROOT :-) *)

  let div ?(c1=1.) ?(c2=1.) h1 h2 =
      if not (fit h1 h2) then fail "mul: different axes";
      let info = make_info "/" (get_info h1) (get_info h2) in
      let h' = dup ~errors:(is_errors h1 && is_errors h2) info h1 in
      let data = get_data h' and 
      data1 = get_data h1 and data2 = get_data h2 in
      let d1 = c1**2. and d2 = c2**2. in (* -- *)
      Array.iterk (fun p data1_p ->
          data.(p) <- (c1 *. data1_p) /. (c2 *. data2.(p));
          let b22 = data2.(p)**2. *. d2 in
          may_errors h1 (fun e1 ->
              may_errors h2 (fun e2 ->
                  may_errors h' (fun e' ->
                      (* -- *)
                      e'.(p) <- d1*.d2*. 
                                ( e1.(p) *. data2.(p)**2. +.
                                  e2.(p) *. data1_p**2. ) /.
                                b22**2.
                 (* d1*d2*(e1*e1*b2*b2 + e2*e2*b1*b1)/(b22*b22) -- *)
          )))
      ) data1;
      h'

  let significance ~sign ~bg =
      if not (fit sign bg) then fail "mul: different axes";
      let info_s = get_info sign and info_b = get_info bg in
      let err_sign = get_errors sign and err_bg = get_errors bg in
      let comment = 
          (if err_sign = None || err_bg = None then "e=sqrt(bin);" else "") in
      let info = { name="sgf("^info_s.name^","^info_b.name^")"; 
          title = comment ^ "significance{ sg = " ^ info_s.title ^ 
                "; bg = "^info_b.title^" }"}
      in
      let data_sign = get_data sign (*and data_bg = get_data bg*) in

      let h = dup ~errors:false info sign in
      let data = get_data h in

      let sq_error h p = ( match get_errors h with
          | Some a -> a.(p)
          | None -> sqrt (get_data h).(p)
      ) in

      Array.iterk (fun p data_sign_p ->
          data.(p) <- data_sign_p /. (sqrt 
                  ((sq_error sign p) +. (sq_error bg p)));
      ) data_sign;
      h


end

module H1 = struct
  open Axis.Bin
  type t = {arr: float array; xa: Axis.t;
        mutable entries: int64; 
        mutable overflow: float; mutable underflow: float; 
        mutable x_nans: float;
        info: info; y_title: string;
        mutable sq_errors: float array option;
    }

  let x_of_bin h k = Axis.c_of_bin h.xa k
  let x1_of_bin h k = Axis.l_of_bin h.xa k
  let x2_of_bin h k = Axis.r_of_bin h.xa k
  let bin_of_x h x = Axis.bin_of_c h.xa x
 
  let x_of_max h = Axis.c_of_bin h.xa (Glob.place_of_max h.arr)
 
  let fill x w h =
     if h.xa &: x
     then (
         let n_bin = bin_of_x h x in
         h.arr.(n_bin) <- h.arr.(n_bin) +. w
     )
     else if x |> h.xa then 
         h.overflow <- h.overflow +. w
     else if x <| h.xa then
         h.underflow <- h.underflow +. w
     else  h.x_nans <- h.x_nans +. w; (* nan !? *)
     h.entries <- Int64.add 1L h.entries;
     ( match h.sq_errors with None -> () | Some _ -> failwith "fill: err" )
  ;;

  let fill_e x w h =
     let err = ( match h.sq_errors with Some e -> e 
       | None -> failwith "fill: no err" ) in
     if h.xa &: x
     then (
         let n_bin = bin_of_x h x in
         h.arr.(n_bin) <- h.arr.(n_bin) +. w;
         err.(n_bin) <- sqrt (err.(n_bin)**2. +. w**2.)
     )
     else if x |> h.xa then 
         h.overflow <- h.overflow +. w
     else if x <| h.xa then
         h.underflow <- h.underflow +. w
     else  h.x_nans <- h.x_nans +. w; (* nan !? *)
     h.entries <- Int64.add 1L h.entries
  ;;
(*  let fill' (x,w) h = fill x w h*)

  let rebin k h =
      Printf.printf "rebinning by %i\n" k; flush stdout;
      let n = h.xa.Axis.n_bins in
      let new_n = n / k in
      let last_bin = new_n * k - 1 in
      let new_axis = Axis.make {n=new_n; x1 = h.xa.Axis.min; 
        x2 = Axis.r_of_bin h.xa last_bin } h.xa.Axis.title in
      let new_arr = Array.create new_n 0. in

      let calc_errors, old_errors, new_errors = ( match h.sq_errors with
        | None -> false, Array.create 0 0., Array.create 0 0.
        | Some a -> true, a, Array.create new_n 0.) in

      for i = 0 to new_n - 1 do
          for j = 0 to k - 1 do
              let old_i = k*i+j in
              new_arr.(i) <- new_arr.(i) +. h.arr.(old_i);
              if calc_errors then (
                  new_errors.(i) <- 
                      new_errors.(i) +. old_errors.(old_i)
              )
          done 
      done;
      let over = ref 0. in
      for l = last_bin to n - 1 do over := !over +. h.arr.(l) done;
      {h with arr=new_arr; xa = new_axis;
      overflow = h.overflow +. !over; 
      sq_errors = if calc_errors then Some new_errors else None}

  exception Found of int
  let rebin_auto h =
      let n = h.xa.Axis.n_bins in
      let max_div = int_of_float (sqrt (float (n+1))) in
      let min_div = try (
          for k = 2 to max_div do 
              if n mod k = 0 then raise (Found k) 
          done;
          n
      ) with Found n -> n in
      rebin min_div h
  ;;


 let h ?(with_errors=false) ?(xname="") ?(yname="") ~(xa: geom) info =
     let cr() =Array.create xa.n 0.0 in
    { arr=cr(); xa = Axis.make xa xname; info=info;
      entries = 0L; overflow = 0.; underflow = 0.; x_nans = 0.;
      y_title=yname; 
      sq_errors=if with_errors then Some (cr()) else None; }

 let h_from_axis ?(with_errors=false) ?(yname="") xa info =
     let cr () = Array.create xa.Axis.n_bins 0.0 in
    { arr=cr(); xa = xa; info=info;
      entries = 0L; overflow = 0.; underflow = 0.; x_nans = 0.;
      y_title=yname; 
      sq_errors=if with_errors then Some (cr()) else None; }

 let set_sq_errors h arr =
     h.sq_errors <- Some arr

 let clear h =
     h.entries <- 0L;
     h.overflow <- 0.;
     h.underflow <- 0.;
     h.x_nans <- 0.;
     Array.iteri (fun i _ -> h.arr.(i) <- 0.) h.arr

 let h_of_array ?(xname="") ?(yname="") arr info =
     let l = Array.length arr in
     let hist = h ~xname ~yname ~xa:{n=l;x1=0.5;x2=0.5+.float l} info in
     for i=0 to l-1 do
         fill (float (i+1)) arr.(i) hist
     done;
     hist

 let analog ?(n_min=1) ?(factor=1) h1 =
     let n = max n_min (factor * (Array.length h1.arr)) in
     h ~xname:(Axis.name h1.xa) ~yname:h1.y_title 
     ~xa:{(Axis.to_geom h1.xa) with n=n} h1.info
 ;;

 include Make_bin_op (struct
     type hist = t let get_info h = h.info
     let dup ?errors info hist = 
         h_from_axis ?with_errors:errors hist.xa info
     let fit h1 h2 = h1.xa = h2.xa
     let get_data h = h.arr let get_errors h = h.sq_errors
     module Array = struct type t = float array type place_info = int
        include Array let iterk = iteri end
 end)

(*
  let divide h1 h2 =
      if h1.xa <> h2.xa then (
          raise (Incorrect_pair "h1: divide")
      ) else (
          let txt = h1.info.name^"/"^h2.info.name in
          let h = h_from_axis h1.xa 
                              {name=txt;title=txt} in
          let n = h.xa.Axis.n_bins in
          ( match h1.y_errors, h2.y_errors with
            | Some e1, Some e2 -> 
                    h.y_errors <- Some (Array.create n 0.)
            | _ -> ());
          for i = 0 to n - 1 do
              h.arr.(i) <- h1.arr.(i) /. h2.arr.(i);
              ( match h.y_errors, h1.y_errors, h2.y_errors with
                | Some e, Some e1, Some e2 ->
                        e.(i) <- sqrt (e1.(i)**2. +. e2.(i)**2.)
                | _ -> ())
          done;
          h
      )
*)
 let x_dots h = Axis.dots h.xa
 let x_left_dots h = Axis.left_dots h.xa
 let range h = Axis.range h.xa

(* let x_dots_last_x h =
     let l=Array.length h.arr in
     let last_x = x_of_bin h l in
     (x_dots h,last_x)
*)
 let integral h =
     Array.fold_left (
         fun acc c -> 
             if not_nan c then acc+.c else acc
     ) 0. h.arr

 let weights h =
     (integral h) +. h.overflow +. h.underflow +. h.x_nans 

 let sum_out h=
     h.underflow +. h.overflow

 let mean h=
     let sumx=ref 0. in
     Array.iteri (fun i content ->
         if not_nan content then
             sumx := !sumx +. content *. (x_of_bin h i)
     ) h.arr;
     let sum = integral h in
     (!sumx /. sum)

 let max_mean_rms h=
     let m=mean h in
     let sumr=ref 0. in
     let _,max = minmax h.arr in
     Array.iteri (fun i content ->
         if not_nan content then (
             sumr := !sumr +. content *. ((x_of_bin h i) -. m) ** 2.;
         )
     ) h.arr;
     let sum = integral h in
     (max, m, sqrt(!sumr/.sum))
(*     Printf.printf "%f %f %f %f\n" max m !sumr sum;*)

let sub (rmin,rmax) h =
    let l=Array.length h.arr in
    let n_1_bin = max 0 (bin_of_x h rmin) in
    let n_last_bin = min (l-1) (bin_of_x h rmax) in
    let n_bins = n_last_bin - n_1_bin + 1 in
    let new_arr = Array.init n_bins (fun i ->
        let n_old = i + n_1_bin in
        if n_old >= 0 && n_old < l 
	then h.arr.(n_old) 
        else 1.) 
    in

    let new_errors = ( match h.sq_errors with
        None -> None
      | Some err -> Some (
              Array.init n_bins (fun i ->
        let n_old = i + n_1_bin in
        if n_old >= 0 && n_old < l 
	then err.(n_old) 
        else 1.))
    ) in

    let underflow=ref h.underflow and overflow=ref h.overflow
    and x_nans = ref h.x_nans in
    Array.iteri(fun i element ->
        let on_left = i<n_1_bin and on_right = i>n_last_bin in
        if on_left then underflow := !underflow +. element;
        if on_right then overflow := !overflow +. element;
        if on_left || on_right then x_nans := !x_nans +. element;
    ) h.arr;
    { h with arr=new_arr; 
    xa = Axis.make {n=n_bins; 
       x1 = x1_of_bin h n_1_bin; x2 = x2_of_bin h n_last_bin} (Axis.name h.xa);
    underflow = !underflow; overflow = !overflow; x_nans = !x_nans;
    sq_errors=new_errors
    }

end

type projection = X | Y
let letter = function X -> "X" | Y -> "Y"



(** 2d histograms *)

module H2 = struct
  open Axis.Bin
    type t = { 
        arr: float array array; xa: Axis.t; ya: Axis.t;
        mutable entries: int64; 
        mutable out: float array array (* 3x3 *);
        mutable xy_nans: float;
        mutable sq_errors: float array array option;
        info: info; }

 let h ?(with_errors=false) ?(xname="") ?(yname="") ~(xa:geom) ~(ya:geom) info =
     let cr() = Array.make_matrix xa.n ya.n 0.0 in
    {arr=cr(); 
    xa = Axis.make xa xname; 
    ya = Axis.make ya yname; 
    entries=0L; out = Array.make_matrix 3 3 0.; xy_nans = 0.; 
    sq_errors = if with_errors then Some (cr()) else None;
    info=info; }

 let h_from_axes ?(with_errors=false) xa ya info =
     let cr () = Array.make_matrix xa.Axis.n_bins ya.Axis.n_bins 0.0 in
     {arr=cr();
    xa = xa; 
    ya = ya; 
    entries=0L; out = Array.make_matrix 3 3 0.; xy_nans = 0.; 
    sq_errors = if with_errors then Some (cr()) else None;
    info=info; }

 include Make_bin_op (struct
     type hist = t let get_info h = h.info
     let dup ?errors info hist = 
         h_from_axes ?with_errors:errors hist.xa hist.ya info
     let fit h1 h2 = h1.xa = h2.xa && h1.ya = h2.ya
     let get_data h = h.arr let get_errors h = h.sq_errors
     module Array = struct type t = float array array type place_info = int*int
        include Matrix end
 end)

 let get_axis h = function X -> h.xa | Y -> h.ya
 let c_of_bin h proj k = Axis.c_of_bin (get_axis h proj) k 
 let x_of_bin h k = Axis.c_of_bin h.xa k
 let y_of_bin h k = Axis.c_of_bin h.ya k
 let xw_of_bin h kx = Axis.w_of_bin h.xa kx
 let yw_of_bin h ky = Axis.w_of_bin h.ya ky
 let x1_of_bin h k = Axis.l_of_bin h.xa k
 let y1_of_bin h k = Axis.l_of_bin h.ya k
 let x2_of_bin h k = Axis.r_of_bin h.xa k
 let y2_of_bin h k = Axis.r_of_bin h.ya k

 let get_bin proj h i j =
     match proj with
     X ->  x_of_bin h i
   | Y ->  y_of_bin h j
(*
 let bin_of proj h p =
     match proj with
     X -> Axis.bin_of_c h.xa p
   | Y -> Axis.bin_of_c h.ya p
*)

 let bin_of proj h p = Axis.bin_of_c (get_axis h proj) p
 let bin_of_x h p = Axis.bin_of_c h.xa p
 let bin_of_y h p = Axis.bin_of_c h.ya p

 let interval h proj k = Axis.interval (get_axis h proj) k
 let lower h proj = Axis.min (get_axis h proj)
 let upper h proj = Axis.max (get_axis h proj)

 let fill x y w h =
     if h.xa &: x && h.ya &: y then (
         let n_bin_x = bin_of_x h x in
         let n_bin_y = bin_of_y h y in
         h.arr.(n_bin_x).(n_bin_y) <- h.arr.(n_bin_x).(n_bin_y) +. w
     )
     else (
         if not_nan x && not_nan y then (
           let sigx = Axis.sgn h.xa x and sigy = Axis.sgn h.ya y in 
             h.out.(sigx).(sigy) <- h.out.(sigx).(sigy) +. w )
         else
             h.xy_nans <- h.xy_nans +. 1.
         );
     h.entries <- Int64.add 1L h.entries
  ;;

  let fill' (x,y) w h = fill x y w h
  let fill'' (x,y,w) h = fill x y w h

 let clear h =
     h.entries <- 0L; Matrix.make_zero h.out; h.xy_nans <- 0.;
     Matrix.make_zero h.arr

 let h_of_array ?(xname="") ?(yname="") arr info =
   let n_i = Array.length arr in
   let n_j = Array.length arr.(0) in
   let hist = h ~xname ~yname 
     ~xa:{n=n_i;x1=0.5;x2=0.5+.float n_i}
     ~ya:{n=n_j;x1=0.5;x2=0.5+.float n_j} info 
   in
     for i=0 to n_i-1 do
       for j=0 to n_j-1 do
         fill (float (i+1)) (float (j+1)) arr.(i).(j) hist
       done
     done;
     hist

 let integral h = Matrix.sum h.arr

 let mean proj h=
     let sump=ref 0. in
     Array.iteri (fun i array ->
         Array.iteri (fun j content ->
             if not_nan content then
                 sump := !sump +. content *. (get_bin proj h i j)
         ) array
     ) h.arr;
     let sum = integral h in
     (!sump /. sum)

 (** returns mean, (mean x, rms x), (mean y, rms y) *)
 let max_mean_rms h=
     let mx=mean X h 
     and my=mean Y h
     in
     let sumx=ref 0. and sumy=ref 0. in
     let max=ref h.arr.(0).(0) in
     Array.iteri (fun i array ->
         Array.iteri (fun j content ->
             if not_nan content then (
                 sumx := !sumx +. content *. ((get_bin X h i j) -. mx) ** 2.;
                 sumy := !sumy +. content *. ((get_bin Y h i j) -. my) ** 2.;
                 if content > !max
                 then max := content
             )
         ) array
     ) h.arr;
     let sum = integral h in
     (!max, (mx, sqrt(!sumx/.sum)), (my, sqrt(!sumy/.sum)) )
(*
  let divide h1 h2 =
      if h1.xa <> h2.xa || h1.ya <> h2.ya then (
          raise (Incorrect_pair "h1: divide")
      ) else (
          let nx = h1.xa.Axis.n_bins in
          let ny = h1.ya.Axis.n_bins in
          let txt = h1.info.name^"/"^h2.info.name in
          let h = h_from_axes h1.xa h1.ya 
          {name=txt;title=txt} in
          ( match h1.z_errors, h2.z_errors with
            | Some e1, Some e2 -> 
                    h.z_errors <- Some (Matrix.create nx ny 0.)
            | _ -> ());
          for ix = 0 to nx - 1 do
            for iy = 0 to ny - 1 do
              h.arr.(ix).(iy) <- h1.arr.(ix).(iy) /. h2.arr.(ix).(iy);
              ( match h.z_errors, h1.z_errors, h2.z_errors with
                | Some e, Some e1, Some e2 ->
                        e.(ix).(iy) <- sqrt (e1.(ix).(iy)**2. +. e2.(ix).(iy)**2.)
                | _ -> ())
            done;
          done;
          h
      )
*)

let sub ~x:x_range ~y:y_range ?(except_last=false) h =
    let xa = match x_range with 
        Some r -> Axis.sub ~except_last h.xa r | None -> h.xa 
    and ya = match y_range with 
        Some r -> Axis.sub ~except_last h.ya r | None -> h.ya in
    let new_h = h_from_axes xa ya 
        {name=h.info.name^", sub";title=h.info.title^" (sub)"} in
    new_h.entries <- h.entries;
    new_h.xy_nans <- h.xy_nans;
    (* TODO: *)
    new_h.out.(0).(0) <- Matrix.sum h.out;

    Array.iteri (fun i a_i ->
        Array.iteri (fun j a_ij ->
            let x = x_of_bin h i and y = y_of_bin h j in
            fill x y a_ij new_h
        ) a_i
    ) h.arr;
    new_h

let projection p h =
    let a = match p with X -> h.xa | Y -> h.ya in
    let h_proj = H1.h_from_axis a 
        {name=h.info.name^" proj";title=h.info.title^" (proj)"} in
    h_proj.H1.entries <- h.entries;
      let sum i =
	( match p with
              X -> sum h.out.(i)
            | Y -> let m = Matrix.transpose h.out in 
		sum m.(i)
	) in
    h_proj.H1.overflow <- sum 2; 
    h_proj.H1.underflow <- sum 0;
    Array.iteri (fun i a_i ->
        Array.iteri (fun j a_ij ->
            let c = match p with
                X -> x_of_bin h i
              | Y -> y_of_bin h j
            in
            H1.fill c a_ij h_proj
        ) a_i
    ) h.arr;
    h_proj

end

type hist=H1d of H1.t | H2d of H2.t

let name hist=
    match hist with
    H1d h -> h.H1.info.name
  | H2d h -> h.H2.info.name

let h1d ?(xname="") ?(yname="") ~xa info = 
    H1d (H1.h ~xname ~yname ~xa info)
let h2d ?(xname="") ?(yname="") ~xa ~ya info = 
    H2d (H2.h ~xname ~yname ~xa ~ya info)

exception Incorrect_axis of string

let get_axis =
    function
        H1d h -> (function X -> h.H1.xa | _ -> raise (Incorrect_axis "Y"))
      | H2d h -> (H2.get_axis h)

let parse_projection = function
     "X" | "x" -> X 
   | "Y" | "y" -> Y
   | s -> raise (Incorrect_axis 
        ("parse_projection: unknown projection '"^s^"'"))

exception Invalid_dimension_hist of string

let only_for_1d txt f = 
    function
     | H1d h -> f h
     | H2d h -> raise (Invalid_dimension_hist (
         "function "^txt^" is available only for 1d histograms"))

let only_for_2d txt f =
    function
     | H2d h -> f h
     | H1d h ->  raise (Invalid_dimension_hist (
         "function "^txt^" is available only for 2d histograms"))
;;
let fill1 x w = only_for_1d "fill1" (H1.fill x w)
let fill2 x y w = only_for_2d "fill2" (H2.fill x y w)

let projection = 
    only_for_2d "projection" (fun h p -> H1d (H2.projection p h))

(*
let fill1 x w hist=
    match hist with
    H1d h -> 
        H1.fill x w h
  | H2d h -> 
        Printf.printf "invalid use of hist 1d \"%s\"\n"
        h.H2.info.title; exit 20

let fill2 x y w hist=
    match hist with
    H2d h -> 
        H2.fill x y w h
  | H1d h ->
        Printf.printf "invalid use of hist 2d \"%s\"\n"
        h.H1.info.title; exit 20
*)

let clear = function
    H1d h -> H1.clear h
  | H2d h -> H2.clear h

let get_1d = function
    H1d h -> h
  | H2d _ -> raise Not_found
;;
let get_2d = function
    H2d h -> h
  | H1d _ -> raise Not_found
;;
let is_1d = function H1d _ -> true | _ -> false
let is_2d = function H2d _ -> true | _ -> false
let select_1d hlist = List.filter is_1d hlist
let select_2d hlist = List.filter is_2d hlist

let init_hist h f=
    let l=Array.length h.H1.arr in
    for k=0 to l-1 do
        h.H1.arr.(k) <- f (H1.x_of_bin h k)
    done;
    h

let init_hist2 h f=
    let lx=Array.length h.H2.arr in
    for kx=0 to lx-1 do
        let ly=Array.length h.H2.arr.(kx) in
        for ky=0 to ly-1 do
            h.H2.arr.(kx).(ky) <- f (H2.x_of_bin h kx) (H2.y_of_bin h ky)
        done
    done;
    h

let dummy_hist () =
    let h=H1.h ~xa:{n=10;x1=0.;x2=10.} {title="Стартовая";name="h"} in
    H1d (init_hist h (fun x -> 1.+. (x-.5.)**2. /. 25.))
let dummy_hist2 () =
    let h=H2.h ~xa:{n=10;x1=0.;x2=10.} ~ya:{n=10;x1=0.;x2=10.} 
            {title="Startted";name="h"} in
    H2d (init_hist2 h (fun x y -> (x/.10.-.0.5)**2. +. (y/.10.-.0.5)**2.))

(** Drawing parameters *)

type log_projection = Log_X | Log_Y | Log_Z
type plot_opt={logscale: log_projection list; grid: projection list}
let is_grid plot_opt proj = List.mem proj plot_opt.grid
let toggle_log axis opt =
    {opt with logscale=
        if List.mem axis opt.logscale 
        then List.filter (fun a -> a <> axis ) opt.logscale
        else axis::opt.logscale
    }

let toggle_grid proj opt =
    {opt with grid=
        if List.mem proj opt.grid
        then List.filter ( (<>) proj ) opt.grid
        else proj::opt.grid
    }


type plot2d=Box|Contour|Grayscale|Color
type plot1d=Lines|Polyline|Columns|Err_y|Line_err
type show_opt={plot1d: plot1d; plot2d: plot2d; color: int }
let standard_plot = {logscale=[];grid=[]}
let standard_show = {plot1d=Lines;plot2d=Box;color=1}
exception Invalid_histogram of string

(*type plot_fun = ?opt:draw_opt -> hist ->  unit*)

let sub_hist ~x ?(y=None) = function
    H1d h ->
       (match x with
         Some r -> 
           H1d (H1.sub r h)
       | None -> 
               raise (Invalid_histogram "1d histogram cannot be divided in y"))
  | H2d h -> 
          H2d (H2.sub ~x ~y h)
;;
let rebin ?(x=None) ?(y=None) = ( function
    H1d h -> 
        H1d ( match x with 
           | Some n -> H1.rebin n h
           | None -> H1.rebin_auto h)
  | H2d h ->
        failwith "rebin h2: not yet implemented"
)

let combine_bin_funs ~f1 ~f2 h1 h2 = 
 (match h1, h2 with
     H1d h1', H1d h2' -> H1d (f1 h1' h2')
   | H2d h1', H2d h2' -> H2d (f2 h1' h2')
   | _ -> raise (Incorrect_pair "cannot combine 1d and 2d hist.")
 )
let div ?c1 ?c2 = 
    combine_bin_funs ~f1:(H1.div ?c1 ?c2) ~f2:(H2.div ?c1 ?c2)
let plus ?c1 ?c2 = 
    combine_bin_funs ~f1:(H1.plus ?c1 ?c2) ~f2:(H2.plus ?c1 ?c2)
let mul = combine_bin_funs ~f1:H1.mul ~f2:H2.mul
let mul_const h c =
  (match h with
     H1d h' -> H1d (H1.mul_const h' c)
   | H2d h' -> H2d (H2.mul_const h' c)
  )

let significance ~sign ~bg =
  (match sign, bg with
    | H1d h1, H1d h2 -> H1d (H1.significance ~sign:h1 ~bg:h2)
    | H2d h1, H2d h2 -> H2d (H2.significance ~sign:h1 ~bg:h2)
    | _ -> raise (Incorrect_pair "cannot combine 1d and 2d hist.")
  )


(** Useful functions *)
module Helpers =
struct
    let mkit t = { name=t; title=t }
end
